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^vj ' Abstract. 

Cauchy-Charactcristic Matching (CCM), the combination of a central 3 + 1 
Cauchy code with an exterior characteristic code connected across a time- 
like interface, is a promising technique for the generation and extraction of 
gravitational waves. While it provides a tool for the exact specification of 
boundary conditions for the Cauchy evolution, it also allows to follow gravitational 
fO ' radiation all the way to infinity, where it is unambiguously defined. 

We present a new fourth order accurate finite difference CCM scheme for 
^Sl 1 a first order reduction of the wave equation around a Schwarzschild black 

^ ■ hole in axisymmetry. The matching at the interface between the Cauchy and 

^i^ ' the characteristic regions is done by transfering appropriate characteristic/null 

fy\ ' variables. Numerical experiments indicate that the algorithm is fourth order 

^~s ^ convergent. As an application we reproduce the expected late-time tail decay for 



o 
o 

u 



the scalar field. 

PACS numbers: 02.60.-x, 02.70.-c, 04.20.-q, 0425. Dm 

1. Introduction 



In the past year there has been remarkable progress in the field of numerical relativity, 
particularly in the simulation of binary black hole systems. Some of today's codes are 
K/i \ capable of tracking the holes for many orbits before the merger and ringdown ^ |2 El ■ 

'V^ ' Since a primary goal of numerical relativity is the production of gravitational wave 

C^ ' templates associated with the merger of black holes, it is important to simultaneously 

make progress with the tools which allow for their accurate generation and extraction. 
The schemes which have been able to push forward the binary black hole 
problem either use spatial compactification ^ , which brings spatial infinity at a finite 
coordinate location, or mesh refinement 01^], which allows the boundary to be placed 
further out by using lower resolution far away from the strong field region. A limitation 
of both methods is that the number of grid-points per wavelength gradually diminishes 
as the radiation propagates away from the sources and therefore the quality of the 
gravitational radiation progressively deteriorates, forcing the extraction to occur at 
a relatively small distance from the final holc:|:. The undesirable effect of a time- like 

I Another reason for extracting radiation at small distances is to avoid boundary effects. 
Nevertheless, the comparison of the real part of r'04 extracted at r = 20M and r = 40M done 
by Baker et al. |^ seems to indicate that the extraction at small distances leads to sufficiently 
accurate results. 
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boundary with a chronological future that intersects the region of interest (where 
physical quantities are calculated), was highlighted in ^|Sj. 

To improve the quality of the outgoing radiation one should aim to keep the 
coordinate distance that the waves have to travel to a minimum, without making their 
speed go to zero. This can be achieved, for example, by employing asymptotically null 
slices El El Ql , or by supplementing the Cauchy evolution with a characteristic 
code [niII3iniIIlIISlIEliniIIHlIIlEniEIlE21ESlElE3Eni, which extends all the 
way to infinity. See Figures Q El and El Both methods resolve incoming waves poorly 
at very large distances, but these are not expected to be significant. The reason for not 
using a single characteristic evolution scheme, without a Cauchy interior covering the 
strong field region, is that such a scheme tends to develop caustics, where neighbouring 
characteristics focus and the coordinate system breaks down. 

In this work we describe a fourth order accurate algorithm for the axisymmetric 
massless Klein-Gordon scalar field around a Schwarzschild black hole, which employs 
a Cauchy evolution scheme, based on space-like slices, and a characteristic scheme, 
based on null slices. The outer boundary of the interior Cauchy region is a time- 
like surface and coincides with the inner boundary of the characteristic region. 
Communication between the two regions is done taking into account the propagation 
of the characteristic variables. At the discrete level the interface is implemented with 
touching grids. We find that when the slices are not simultaneous, this particular 
multi-grid approach is better suited than overlapping grids, as it avoids complications 
due to causality. The radial coordinate of the characteristic region is compactified in 
order to reach future null infinity. 

The paper is organised as follows. In SectionElwe describe the continuum scalar 
field initial-boundary value problem in the Cauchy and characteristic regions and 
discuss the interface treatment. In SectionElwe provide details of the discretization 
method and in Section El we describe the tests carried out to validate the code, 
including the simulation of the late-time tail decay. 

2. Klein-Gordon scalar field around a Schwarzschild black hole 

In this section we describe the continuum problem, giving the explicit form of the 

equations, and the interface treatment. The Schwarzschild line element in Kerr-Schild 

coordinates takes the form 

2M 

ds^ = -dt^ + dr^ + (dt + drf + r'^dn^, (1) 

r 

where M is the mass of the black hole and dil^ — d9^ + sin 9^d(fP . The surface defined 

by 

t = r + AM ln(r - 2M) + const (2) 

is null and its normal, provided it belongs to the future null cone, points toward 
increasing values of r. It is convenient to introduce the shorthand r = r±2Af . Using 
\1/ — r<f>, in the Cauchy region the wave equation on the Schwarzschild background, 
V^V^$ = 0, takes the form 

In Bondi coordinates the Schwarzschild line element is 
ds^ = -IH dfi - 2 du dr 



2M de{sin0dg-^) 
r+r^ rr+ sin 


(3) 


ment is 




7-'dn^ 


(4) 
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Figure 1. Most simulations in numerical relativity are carried out on a domain 
bounded by two space-like slices, So a-nd Si (the initial and final time), and a time- 
like artificial boundary T. Since it is not known how to specify local boundary 
conditions which do not introduce spurious radiation or reflect outgoing waves, 
the accuracy of the solution is limited to the future domain of dependence of So 
(green region) . The solution in the grey region, which is causally connected to the 
boundary, is inaccurate. In this situation one is forced to extract gravitational 
radiation relatively close to the sources. 



and it is related to (^ via the coordinate transformation 
u = t-r-AMln{r-2M), 



which imphes 



dt = du, 



(5) 
(6) 
(7) 
(8) 

(9) 
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Figure 2. When asymptotically null slices are used, as in 1271 . the outer boundary 
T can be placed far away from the sources and is approximately null (the incoming 
modes travel very slowly). Gravitational radiation can be accurately extracted 
at much larger distances and with lower computational cost than with standard 
space-like slices. See Table I of 1271 for a comparison of computational costs and 
errors with different slices. 






2M 



de = dg, 
d^ = dr. 



(10) 

(11) 
(12) 



Where there is no cause for confusion we will not distinguish between {f, 0, (f)} and 

{r,e,^}. 

The Klein-Gordon equation in the characteristic region becomes 



2r ^ , r- ^n 2M ^ 2M ^^-(sine'as*) , , 

= -—dudf'^ + — 9|^ + — -9f ^ - — ^^ + ' , . ' ■ 13 
J.+ j,+ j,j.+ j,+j,^ rr+ sint/ 
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Figure 3. By combining a Caucliy and a characteristic scheme one solves for 
the entire spacetime of interest, including J^+ (future null infinity), where the 
radiation is unambiguously defined. Initial data is specified at So and A^o and the 
two schemes are connected by a time-like interface /. Although in principle one 
could use the asymptotically null approach up to null infinity, experiments with 
the wave equation in spherical symmetry suggest that the scheme is less effective 
than that which uses a time-like boundary at very large distances 1271 . 



Note the lack of second partial derivatives 9,^, reflecting the fact that the u = const 
surface is characteristic. 



2.1. Axisymmetry and first order reductions 

We focus on axisymmetric solutions of the equations Q and l(T^ , as we do not expect 
additional difficulties to arise in the general 3-dimensional case. This implies 9^$ = 
and regularity demands that, across the axis 9 = and = tt, the scalar field is an 
even function of 9. 
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Since the discretization of second order systems tends to be more subtle than that 
of fully first order systems |2H1 HSl QOl 1^ 112] j we perform a reduction to first order 
and discretise the resulting system. In the Cauchy region we introduce the auxiliary 
variables T = dt'^, R = dr'i' and G = de'^- The reduced system is 

4A/^ 7'-^ 2 A/ 2M 2M ^e (sin 6*9) , 

J.+ j.-t- j,^+ j,^+ j.zj,+ rr+ sine/ 

dtR = drT, (15) 

5*6 = deT, (16) 

dt-^ = T. (17) 

We also have the constraint C = drQ — dgR ~ 0, which propagates with zero speed 
and therefore does not introduce complications at the boundary. This system is 
symmetrizable hyperbolic provided that r > 2M |33| . It has non trivial characteristic 
speeds 1 and —r^/r^ in the radial direction and the corresponding characteristic 
variables are 

r+ 
w'^ = —T±R. (18) 

r 

To reduce equation H13|) to first order we introduce R = d,-'^, V = 2r/7'^9„^— 9^ ^ 
and O = dg'i'. The reason for this choice of auxiliary variables is that df and 
2r/r~du — df are null vectors. We will refer to R and V as null variables. The 
reduced system is 

5.y =-^^-4^* + ^!!!^, (19) 

rr r^r rr s\nO 

^ ~ r- ^ ~ 2M ~ 2M d^ismOQ) , , 

2duR^—dfR+^R^^^^+ '\ . ' , (20) 

r r^ r'^ r^ sm 

2d^Q^'--{d~,V + d^R), (21) 

2du^ = —{V + R). (22) 

r 

To reach ^"'" we compactify the characteristic radial coordinate r by introducing 
a new independent variable ^ such that 

f = n + -{i^-£.^)t^n(^^^\ (23) 

In the last equation r = ri and ^ = ^i represent the location of the interface boundary 
in the Cauchy and characteristic coordinate systems, respectively. The surface ^ = ^oo 
corresponds to future null infinity. Effectively, this amounts to replacing in system 
Cnj"|23 the partial derivatives with respect to f with 



* = ^a.=-^^-^=-e|ja, (24) 

There are coefficients in equation (|19|l that contain indetermined forms, which can be 
explicitly evaluated using the limit 

hm f2p = 4(?--^^)'- (25) 

?^«oo or -K^ 

'SNc typically use r^ = (t = lOM and ^oo = 20i\f . 
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2.2. Interface treatment 

A necessary condition for well-posedness for the interior Cauchy problem is that at 
the boundary and at the interface data is prescribed to the incoming characteristic 
variables. If no data is specified to these variables, then one can expect to lose 
uniqueness. The characteristic region can provide accurate data for the incoming 
mode by setting 

W+\r=r,^Rk=i,. (26) 

On the other hand, the characteristic problem requires data for the outgoing null 
variable V. This is done by setting 

V\^=^^^[T~R]r^r,. (27) 

2.3. Treatment of future null infinity 

Future null infinity corresponds to the boundary surface ^ = ^oo. Here, as a result of 
the compactification, the speed of the incoming variable R approaches zero, so there 
are no incoming modes (the surface is null). No data is required at this boundary and 
its numerical treatment is particularly simple. See subsection 13.51 

3. Finite difference scheme 

Stable finite difference approximations of the initial-boundary value problem for the 
wave equation in the Cauchy region (not coupled to the characteristic problem) are 
available. For example, by approximating the spatial derivatives using finite difference 
operators that satisfy the summation by parts rule [3^, it is possible to construct 
schemes which conserve a discrete energy [6'6\ I35| . a positive definite quadratic form 
of the dependent variables. Boundary conditions can then be imposed cither using 
Olsson's method [HniEZIEHl or using penalty terms P^HTIET] . 

It is possible to obtain energy estimates in the characteristic region |42l 1431 E| , 
and in turn for the coupled Cauchy-characteristic problem. Ideally, one would like to 
achieve similar estimates at the discrete level. Although we have been able to formally 
extend these estimates to the discrete case, the semi-discrete problem that we consider 
fails to admit solutions for generic initial data. 

Formally, one can obtain discrete estimates for the semi-discrete system 

A^duw + A'D^w + A^DyW + A^D^w == 0, (28) 

where w^ which is a shorthand for Wijk(u), is a vector- valued grid- function having 
the dynamical variables as components and the variables x, y and z are discretised 
with indices i, j and k. Furthermore, u is left continuous, y and z arc periodic, the 
matrices are constant and symmetric and A" is singular. The operators I?^;, Dy and 
Dz are discrete finite difference operators. The operators Dy and D^ are centered 
finite difference approximations of dy and dz . The operator D^ approximates dx and 
satisfies the summation by parts rule§ 

N N 

'^{DxWi)^Viaih + ^ Wi{DxVi)a^h = wnVn - wqVq, 

?=0 i=0 

§ The summation by parts rule is the discrete analogue of the integration by parts rule: J w'{x)v(x)+ 
w{x)v'{x) dx = w{b)v{b) — w{a)v{a). 
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where Ci are appropriate weights (cxphcit examples of weights and operators can be 
found in [SI]) and h is the spatial mesh spacing. System (|28|l represents the semi- 
discrete approximation of H19|l - H22(l . Without loss of generality we can assume that 
A" is diagonal and therefore at least one equation in system (|28ll does not contain 
derivatives with respect to u. Multiplying on the left by w'^ and using the fact that 
the matrices are constant and symmetric, we get 

{duwfA'^w + w'^A^'duW + {D^wfA'^w + w'^A'D^w 

+ {DyW^A^w + uFAyOyW + [D^wfA'^w + uFA'^D^w = 0. (29) 

Summing over the indices corresponding to x, y and z and integrating in u gives the 
equality 



■ N 

^ wfA^'w.a.h 



[wfA-w,]]Jdu = 0, (30) 

Li=0 

where [f{u)]'^Zul = /(^i) ~ /(""o), [fi\l=o = In - /o, the index i corresponds to 
the X direction, and the indices corresponding to y and z are suppressed. Despite the 
estimate, a solution of (|28() does not exist. This can be seen by looking at the equation 
that does not involve derivatives in the u direction. The system DxW ~ F, for generic 
F, is overdetermined because the operator Dx is singular. 

Since the summation by parts approach does not appear to work, our guiding 
principle will be to discretise the problem in a way consistent with the continuum 
problem[|, for which estimates do exist. Numerical experiments will be crucial for 
establishing convergence of the resulting scheme. 

3.1. Cauchy grid 

The discretization of the Cauchy problem is done in a standard way: the derivatives 
are approximated with centered fourth order accurate finite difference operators 

D^ U - \d+d\ , (31) 

where Dqw, == {w^+l-w^-l)/{2h), D+Wi = {wi+l~w^)/h, and ZJ^mj = {wi-w.i^i)/h. 
When applied to a grid- function Wi, operator H31|l gives 

{-Wi+2 + Bwi+i - Swj-i + mi_2)/(12/i). 

The boundaries, the interface and the axis of symmetry are implemented using 
ghost points. To eliminate the black hole singularity from the domain we use 
excision. At the boundary the ghost points are populated^ using extrapolation and by 
overwriting the incoming characteristic fields, if any are present, at the boundary point. 
See Appendix B.2 of |^ and 0^. Near the axis of symmetry these are populated 
using the regularity conditions, i.e., using the fact that ^ is an even function of 9 
across the axis. The interface boundary is treated as if it was an outer boundary, with 
the only exception that the incoming mode is provided by the characteristic evolution. 

We also add some artificial dissipation to the scheme of the form a-h^{D^D^)^, 
explicitly 

a{wt+3 - 6wi+2 + I5wi+i - 20wi + 15wi„i - 6wi^2 + Wj-a)/^, 

II A finite diflference scheme is said to be consistent if solutions to the continuum problem satisfy the 

scheme to an order greater or equal to 1 in the mesh spacing. 

^ Populating ghost points means assigning the values of the dynamical variables at these points. 



Exact boundary conditions in numerical relativity ... 9 

with a — 0.075 [461 and use fourth order Runge-Kutta to integrate in time with 
a Courant factor of 0.5. Although we are investigating a hnear problem, the 
partial differential equations have variable coefficients. A certain amount of artificial 
dissipation is desirable for stability and will most likely be necessary in the non-linear 
case. 

3.2. Characteristic grid 

The discretization of the evolution part of the system, Eqs. H20|l - H22() . is done as in 
the Cauchy case. Equation H19|l treated as an ordinary differential equation in f, the 
last two terms in the right hand side taking the role of forcing terms. It is advanced 
in the r direction by using fourth order Rungc-Kutta (with a "time" step equal to 
the radial mesh spacing). This requires the evaluation of the coefficients and forcing 
terms of H19(l at intermediate time steps. The missing data is generated with fourth 
order interpolation 

U^+l/2 = i-U,-i + 9U, + 9C/,+i - C/,+2)/16. (32) 

Note that to start the integration we need the fields at the grid-point i = — 1 . We use 
fourth order extrapolation to populate this point. 

The integration of V along the f direction is done at each intermediate time step 
of the global time integrator. Artificial dissipation is added to the evolution in the u 
direction. 

3.3. Axis of symmetry 

We use ghost points to treat the axis 6* = 0, tt. The dynamical variables ^P, T, R, 
R and V are even functions of 9 along the axis, and Q and Q are odd functions. In 
particular, we impose that = G = on the axis. 

3.4. Interface treatment 

Data is communicated at the interface at each intermediate time step: the Cauchy 
code requires the incoming variable w'^ , Eq. ]2<j\i , and the characteristic code needs the 
outgoing null variable V, Eq. H27|l . The Cauchy and characteristic grids are touching, 
rather than overlapping. This avoids complications with causality which would arise 
if the u = const, slice intersected the domain of dependence of the t = const slice. 

At the interface boundary we implement condition l|26(l by solving for T and R 
the system of equations 

(33) 

—T -R= —Told - Raid. (34) 

r r 

3.5. Treatment of future null infinity 

The only difficulty at the outer boundary of the characteristic grid is that some 
coefficients contain indetermined forms "0 ■ 00". We use the limit H25|l to evaluate 
these. We populate the ghost points by extrapolating all fields. As there are no 
incoming modes no further action is taken. 



T + R = 
r 


-R, 


T+ 


r+ 
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3.6. The algorithm 

In this subsection we list the basic steps of our algorithm. 

(i) Give initial data to all fields in the Cauchy grid (^ and T can be given arbitrarily, 
but R and have to obey their definition constraints). Give initial data to all 
fields on the characteristic grid (^ can be given arbitrarily, R and 8 have to 
satisfy their definition constraints), with the exception of the null variable V. 

(it) Impose regularity on the axis, impose interface boundary conditions 126|) and 127|) . 
Integrate Eq. (|19|) using fourth order Runge-Kutta (use fourth order Lagrange 
interpolation to populate the coefficients and forcing terms of the differential 
equation). Extrapolate all fields in the radial direction to populate the ghost 
points. Impose regularity again. 

(iii) Evaluate the right hand sides of the evolution equations (because we use ghost 
points this can be done at every interior, and boundary/axis/interface point) and 
calculate the solution at the next intermediate time step. 

(iv) Repeat from point ^ until a full fourth order Runge-Kutta time step is taken. 

(v) Repeat from point ((nj and take as many Runge-Kutta time steps are needed to 
reach the final time. 

4. Numerical tests 

4.1. Convergence test 

To test for convergence we set M = 0, construct an exact solution by translating the 
spherically symmetric solution $ = /(t — r)/r along the axis by an amount i5, 

* = ?■$ = ^/(t-f), (35) 

r 

where f^ = r^ — 2r6cos9 + S^. The expressions for T, R, Q and R, V, Q can be 
readily computed. We monitor how the L2-iiorm of the error between the numerical 
and exact solution, 

1/2 

\W - M(cxact)ll = X]^^J ~ "(exact) )^^r/le 

scales with resolution. Since the exact solution is known, we only need to use two 
resolutions (coarse and fine) to test for convergence. In the expression above hr and 
kg represent the mesh spacing and the sum is extended to all grid-points of the Cauchy 
and charateristic grids. 

The coarse resolution is 200 x 200 in each grid. The fine resolution is obtained 
by doubling the number of grid-points in each direction. The Cauchy slices extend 
from r = 2, where an artificial boundary is introduced to avoid the r = coordinate 
singularity, to r = 10, the location of the spherical interface boundary. We use a time 
step At = Au ~ 0.02 for the coarse case and half it for the fine case. The test confirms 
that the scheme is indeed fourth order convergent. See figure 01 
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Figure 4. This plot illustrates the experimental convergence rate of the scheme 
as a function of time. The initial pulse, which is of compact support, completely 
leaves the domain shortly before t = 15. The numerical errors left in the 
computational domain converge away at the correct rate. The inset shows the 
L2-norms of the errors at the different resolutions. 



4-2. Tail decay 

Having established that the scheme is fourth order accurate we want to determine 
whether it can accurately produce the late time tail decays 071021 0^1 • The field ^ is 
expected to decay as t~^ (along a time-like curve and assuming spherically symmetric 
initial data), because of the scattering of the field off the curvature of spacetime 
asymptotically far from the black hole. We set M = 1, give initial data to 'J consisting 
of a spherically symmetric pulse with compact support, 

cos'^(r-4), |r-4| < tt/2, 
0, otherwise. 



*(0,r,6l) = 



and set its time derivative equal to zero. We use r.i ~ ^^ = 10 and ^oo = 20. We 
monitor the value of ^ at three different locations: r = 6, r = 16 (^ = 15) and £, = ^tx3 
(i.e., at J^"*"). Using 400 grid-points in the radial direction in both grids and very 
few in the angular direction, we obtain good agreement with the expected result. See 
figure O The errors in the angular direction arc of round-off order, i.e., 10~^^ times 
smaller than the solution. 

5. Conclusion 



To study gravitational radiation at future null infinity, where it is unambiguously 
defined, numerical schemes have to evolve the fields up to that surface. An approach 
which would allow for this is Cauchy-characteristic matching (CCM), whereby 
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1000 



Figure 5. Amplitude of the scalar field as measured by observers at fixed radial 
coordinate r = 6M, r = 16M and at ^+. Whereas the first observer is in the 
Cauchy region, the other two are in the characteristic region. Note that the decay 
rate measured by the observer at infinity differs from that of the observers at a 
finite distance. At t = 800M the ratio (i{ln'I')/d(lni) is -2.92, for the r = 6M 
and r = 16M observers, and —1.97, for the observer at ,J^~^ . 



Cauchy and characteristic evolutions are matched across a time-hke world-tube. 
This approach presents many advantages over other boundary conditions, such as 
maximally dissipative or absorbing boundary conditions. CCM makes it possible, for 
example, to compute accurate waveform and polarisation properties at infinity. It 
is a highly computationally efficient technique which can eliminate the undesirable 
reflections and inaccuracies that other boundary conditions would introduce \l'6\ . 

In this work we described a fourth order accurate algorithm for a Cauchy- 
characteristic formulation of the wave equation written in first order form. By 
compactifying the radial coordinate in the characteristic region we can evolve all the 
way to ^~^ . We carried out convergence tests showing that the scheme is fourth order 
convergent and were able to reproduce the expected tail decay rate. The interface was 
implemented using touching grids, which share a common smooth boundary, and by 
communicating the relevant characteristic/null variables. 

The success for the scalar field case is clearly only a necessary condition for the 
fully nonlinear case. The nonlinear extension of this work will require a careful choice 
of coordinate conditions and attention to the propagation of constraints )5()ll51ll^l5Hj . 
It will be crucial to establish how constraint violations propagate across the domain, so 
that the interface treatment can be done consistently. In three dimensions one would 
still require a smooth interface boundary implemented through either overlapping grids 
| 46l 1541 I55j , which work well with simultaneous slices and can also allow for moving 
boundaries, or touching grids ISEIEZIISHI- 
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